;********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;load library for  wetbulb_stull
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/heat_stress.ncl"
;The *area_conserve_remap *and *area_conserve_remap_Wrap* have been *deprecated*.
;The* ESMF regrid* <http://www.ncl.ucar.edu/Applications/ESMF.shtml> library is the recommended replacement. It uses a superior methodology.
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
;********************************************
 begin

;set baseline
 dirPath = "~/source_data/"
 dirPath2 = dirPath+"/outdoor/"  ;outdoor sun/shade/dark scenarios as in Table 1
 dirPath3 = dirPath+"/work_noSwLimit/"; dark* scenario 
 dirPath4 = dirPath+"/rest_noSwLimit/" ;  dark** scenario

 mo  = (/"ACCESS-CM2","AWI-CM-1-1-MR","BCC-CSM2-MR","CMCC-CM2-SR5","EC-Earth3","GFDL-ESM4",\
              "IITM-ESM", "KIOST-ESM","MIROC6","MIROC-ES2L","MPI-ESM1-2-HR","MRI-ESM2-0"/); 12 models
 case = (/"ssp585" /)
 ncase = dimsizes(case);
 nmod = dimsizes(mo)

 syr = 1990; including historical data to identify 1C warming
 fyr = 2099;
 csyr = sprinti("%0.4i",syr);
 cfyr = sprinti("%0.4i",fyr);
 nyr = fyr - syr + 1;
 ndays = nyr*365

;read dimension for 1 deg data
  dimPath=dirPath+"/domain/"
  landfrac = "sftlf_fx_360x180_1deg_gn.nc"
  area = "areacella_fx_360x180_1deg_gn.nc"

  fdim  = addfile(dimPath+landfrac, "r")
  fdim2 = addfile(dimPath+area, "r")
  lat := fdim->lat
  lon := fdim->lon
  nlat = dimsizes(lat);
  nlon = dimsizes(lon);
  garea = fdim2->areacella  ;m^2
  lfrac = fdim->sftlf/100; %
  wgt := garea*1e-6*lfrac ; sqrkm
  wgt@_FillValue = lfrac@_FillValue; make sure FillValue is set correctly for wgt
  copy_VarCoords(garea, wgt)
  printVarSummary(wgt)
  print("global total land area: "+sum(wgt)+" sq km")



  ;----data from Data.Outdoor_Gday_impacts_gridcell-ens-stat-vcbc.ncl
  ;spatial statistics and extent of impact
  csv_file1 = dirPath2+"Gday_hot1pct_sun_warminglevels_10yravg_ens_median-vcbc-outdoor.csv"   ;sun scenario
  csv_file2 = dirPath2+"Gday_hot1pct_diff_warminglevels_10yravg_ens_median-vcbc-outdoor.csv"   ;shade scenario
  csv_file3 = dirPath2+"Gday_hot1pct_zero_warminglevels_10yravg_ens_median-vcbc-outdoor.csv" ;dark scenario

  ;TW based impacts
  csv_file10 = dirPath2+"TWday_hot1pct_warminglevels_10yravg_ens_median-vcbc-outdoor.csv"

  ;dark* scenario in Table 1
  csv_file31 = dirPath3+"Gday_hot1pct_zero_warminglevels_10yravg_ens_median-vcbc-indoor.csv"; 
  ;dark** scenario in Table 1
  csv_file32 = dirPath4+"Gday_hot1pct_zero_warminglevels_10yravg_ens_median-vcbc-indoor.csv"; 


;---Read in file as array of strings so we can parse each line
   lines  := asciiread(csv_file1,-1,"string")
   lines2 := asciiread(csv_file2,-1,"string")
   lines3 := asciiread(csv_file3,-1,"string")
   lines10 := asciiread(csv_file10,-1,"string")
   lines31 := asciiread(csv_file31,-1,"string")
   lines32 := asciiread(csv_file32,-1,"string")
   nlines = dimsizes(lines)-1   ; First line is a header
   delim = ","

  ;---skip the header line
   lines  := lines(1:)
   lines2 := lines2(1:)
   lines3 := lines3(1:)
   lines10 := lines10(1:)
   lines31 := lines31(1:)
   lines32 := lines32(1:)

;---Read fields (columns)
   levs =  tofloat(str_get_field(lines,1,delim))
   ;--select the warming levels to plot
   levels := fspan(1,4,31); 1 to 4 degC by 0.1 degC
   levid = get1Dindex(levs,levels)
   nlev = dimsizes(levels)

   Gday999g_sun_median =   tofloat(str_get_field(lines(levid),2,delim))
   Gday99g_sun_median =   tofloat(str_get_field(lines(levid),3,delim))
   Gday95g_sun_median =   tofloat(str_get_field(lines(levid),4,delim))
   Gday999g_sun_high =   tofloat(str_get_field(lines(levid),5,delim))
   Gday99g_sun_high =   tofloat(str_get_field(lines(levid),6,delim))
   Gday95g_sun_high =   tofloat(str_get_field(lines(levid),7,delim))
   Gday999g_sun_low =   tofloat(str_get_field(lines(levid),8,delim))
   Gday99g_sun_low =   tofloat(str_get_field(lines(levid),9,delim))
   Gday95g_sun_low =   tofloat(str_get_field(lines(levid),10,delim))

   Gday999g_diff_median =   tofloat(str_get_field(lines2(levid),2,delim))
   Gday99g_diff_median =   tofloat(str_get_field(lines2(levid),3,delim))
   Gday95g_diff_median =   tofloat(str_get_field(lines2(levid),4,delim))
   Gday999g_diff_high =   tofloat(str_get_field(lines2(levid),5,delim))
   Gday99g_diff_high =   tofloat(str_get_field(lines2(levid),6,delim))
   Gday95g_diff_high =   tofloat(str_get_field(lines2(levid),7,delim))
   Gday999g_diff_low =   tofloat(str_get_field(lines2(levid),8,delim))
   Gday99g_diff_low =   tofloat(str_get_field(lines2(levid),9,delim))
   Gday95g_diff_low =   tofloat(str_get_field(lines2(levid),10,delim))

   Gday999g_zero_median =   tofloat(str_get_field(lines3(levid),2,delim))
   Gday99g_zero_median =   tofloat(str_get_field(lines3(levid),3,delim))
   Gday95g_zero_median =   tofloat(str_get_field(lines3(levid),4,delim))
   Gday999g_zero_high =   tofloat(str_get_field(lines3(levid),5,delim))
   Gday99g_zero_high =   tofloat(str_get_field(lines3(levid),6,delim))
   Gday95g_zero_high =   tofloat(str_get_field(lines3(levid),7,delim))
   Gday999g_zero_low =   tofloat(str_get_field(lines3(levid),8,delim))
   Gday99g_zero_low =   tofloat(str_get_field(lines3(levid),9,delim))
   Gday95g_zero_low =   tofloat(str_get_field(lines3(levid),10,delim))

   TWday999g_median =   tofloat(str_get_field(lines10(levid),2,delim))
   TWday99g_median =   tofloat(str_get_field(lines10(levid),3,delim))
   TWday95g_median =   tofloat(str_get_field(lines10(levid),4,delim))
   Gtw999g_median =   tofloat(str_get_field(lines10(levid),5,delim))
   Gtw99g_median =   tofloat(str_get_field(lines10(levid),6,delim))
   Gtw95g_median =   tofloat(str_get_field(lines10(levid),7,delim))
   TWday999g_high =   tofloat(str_get_field(lines10(levid),8,delim))
   TWday99g_high =   tofloat(str_get_field(lines10(levid),9,delim))
   TWday95g_high =   tofloat(str_get_field(lines10(levid),10,delim))
   Gtw999g_high =   tofloat(str_get_field(lines10(levid),11,delim))
   Gtw99g_high =   tofloat(str_get_field(lines10(levid),12,delim))
   Gtw95g_high =   tofloat(str_get_field(lines10(levid),13,delim))
   TWday999g_low =   tofloat(str_get_field(lines10(levid),14,delim))
   TWday99g_low =   tofloat(str_get_field(lines10(levid),15,delim))
   TWday95g_low =   tofloat(str_get_field(lines10(levid),16,delim))
   Gtw999g_low =   tofloat(str_get_field(lines10(levid),17,delim))
   Gtw99g_low =   tofloat(str_get_field(lines10(levid),18,delim))
   Gtw95g_low =   tofloat(str_get_field(lines10(levid),19,delim))
   hc999g_med =   tofloat(str_get_field(lines10(levid),20,delim))
   hc99g_med =   tofloat(str_get_field(lines10(levid),21,delim))
   hc95g_med =   tofloat(str_get_field(lines10(levid),22,delim))

   ;--add scenario without sweat limit but with moderate work MH=300/1.7
   Gday999g_zero_median_noSw =  tofloat(str_get_field(lines31(levid),2,delim))
   Gday99g_zero_median_noSw =   tofloat(str_get_field(lines31(levid),3,delim))
   Gday95g_zero_median_noSw =   tofloat(str_get_field(lines31(levid),4,delim))

   ;--add scenario without sweat limit and resting (MH=100/1.7)
   Gday999g_zero_median_noSw_rest =  tofloat(str_get_field(lines32(levid),2,delim))
   Gday99g_zero_median_noSw_rest =   tofloat(str_get_field(lines32(levid),3,delim))
   Gday95g_zero_median_noSw_rest =   tofloat(str_get_field(lines32(levid),4,delim))



  ;--calculate the nonlinear scale to align Gtw on the left Y axis to TW on the right Y 
  cp=1010 ;specific heat capacity of air (Joules/kg/°C)
  pa=101000 ;air pressure (Pa);  
  Ts =35 ;default for TW and Gtw (Eq. 10/14) as used by Sherwood and Huber (2010)
  lambda = latent_heat_water(Ts, (/0, 0/), 1, False)

  print("ensemble median TW_day of hottest 1% grid cells at each warming level: "+TWday99g_median+" at "+levels) 
  print("ensemble median hc_day of hottest 1% grid cells at each warming level: "+hc99g_med+" at "+levels) 
  print("mean hc_day of hottest 1% grid cells across warming levels: "+avg(hc99g_med))

  ;--use dynamic hc and k corresponding to TW of 25-38 to estimate the nonlinear relationship by Eq.14
  ;Eq. 14: Gtw = -hc(1+lambda/cp*delta)(Ts-Tw) = k(Ts-Tw)
  ;use regularly spaced Gtw to estimate Tw: Tw = Ts - Gtw/k

  ;left and right Y axis limits
  Tw_min = 25
  Tw_max = 38

  ;when delta (k) is evaluated at (TW+Ts)/2
  ;Tw_right = fspan(Tw_min,Tw_max, (Tw_max-Tw_min)*10+1); increment by 0.1
  Tw_right = fspan(Tw_min,Tw_max, (Tw_max-Tw_min)*5+1); increment by 0.2
  print(Tw_right)
  Gtw = new(dimsizes(Tw_right), float)
  do i=0,dimsizes(Tw_right)-1
    TW = Tw_right(i)

    if (TW .ge. min(TWday99g_median-1) .and. TW .le. max(TWday99g_median)) then ; *make sure the scale is monotonically increasing!
      tDif = abs(TW - TWday99g_median)
      ind1 = minind(tDif); closest index
      print("TW="+TW+" and Gtw="+Gtw99g_median(ind1)+" at index "+ind1)
      hc = hc99g_med(ind1)
    else; use mean hc outside of TW range
      hc = avg(hc99g_med)
    end if 
    delta_q := (0.622/pa)*1000*(2508*exp(17.3*(TW+Ts)/2 /((TW+Ts)/2 + 237.3))/((TW+Ts)/2 + 237.3)^2)
    k := -hc*(1+lambda*delta_q/cp); use ensembe median hc
    Gtw(i) = k*(Ts-TW)
  end do

  G_scale = Gtw; used as nonlinear scale on left Y axis
  Gtw_min = min(Gtw)
  Gtw_max = max(Gtw)
  print("max-min range of Gtw on the left Y axis corresponding to Tw range of "+Tw_min+"-"+Tw_max+" : "+Gtw_min+"-"+Gtw_max)
  printMinMax(Gtw99g_median,1)
  printMinMax(TWday99g_median,1)
  print(G_scale); *make sure the scale is monotonically increasing!


  ;-----Fig.5ab------

  ; define a polygon centered the ensemble median 
  ;================================================
   xx = levels ;show data from 1 to 4 degC by 0.1 degC

   xp    = new( (/2*nlev/), float )
   yp1   = new( (/2*nlev/), float )
   yp2   = new( (/2*nlev/), float )
   yp3   = new( (/2*nlev/), float )
   yp4   = new( (/2*nlev/), float )

   xp(0:(nlev-1)) = xx
   xp(nlev:(2*nlev-1)) = xx(::-1) ; reverse

   ;ensemble CI for top 1% grid cells
   yp1(0:(nlev-1)) = Gday99g_sun_high ; uppper bound
   yp1(nlev:(2*nlev-1)) = Gday99g_sun_low(::-1) ; lower bound
   yp2(0:(nlev-1)) = Gday99g_diff_high ; uppper bound
   yp2(nlev:(2*nlev-1)) = Gday99g_diff_low(::-1) ; lower bound
   yp3(0:(nlev-1)) = Gday99g_zero_high ; uppper bound
   yp3(nlev:(2*nlev-1)) = Gday99g_zero_low(::-1) ; lower bound
   yp4(0:(nlev-1)) = TWday99g_high ; uppper bound
   yp4(nlev:(2*nlev-1)) = TWday99g_low(::-1) ; lower bound

   print("OK")

 ;---draw line plots
 imagePath = dirPath+"/../image/" 
 image1 = imagePath+"/Fig.3"
 wks1 = gsn_open_wks("pdf",image1)

 res = True
 res@gsnDraw              = False
 res@gsnFrame             = False

 ;res@vpXF                  = 0.1;0.15
 res@vpWidthF             = 0.34; 0.6 ;plot width
 res@vpHeightF            = 0.3; 0.3   ;plot length

 res@tiXAxisFontHeightF   = 0.014;0.022
 res@tiYAxisFontHeightF   = 0.014;0.022
 res@tmXBLabelFontHeightF = 0.012
 res@tmYLLabelFontHeightF = 0.012
 res@tmYRLabelFontHeightF = 0.012
 res@gsnStringFontHeightF = 0.014
 res@tmXTOn               = False
 ;res@tmXBMinorOn          = False
 ;res@tmBorderThicknessF    = 1.5

 res@pmLegendDisplayMode    = "Never";"Always"
 res@pmLegendSide           = "Top"; reference axis
 res@pmLegendParallelPosF   = .38; left-right adjustment
 res@pmLegendOrthogonalPosF = -.31;
 res@pmLegendWidthF         = 0.4;0.1
 res@pmLegendHeightF        = 0.04;0.4
 res@lgLabelFontHeightF     = .015
 res@lgPerimOn              = False
 res@lgOrientation = "horizontal"
 res@lgLabelAngleF = 300

 cmap = read_colormap_file("srip_reanalysis");("radar")
 cmap = cmap(::-1,:) ; reverse color map
 cmap := cmap((/1,2,3,4,5,6,7,8,10,11,12,13,15,16,18,0/),:); first black, last gray for ensemble median
 res@tiYAxisString   = "G~S~day~N~ or G~B~TW~N~~S~day~N~ (W m~S~-2~N~)"
 res@tiXAxisString   = "Warming since pre-industrial (~S~o~N~C)" 
 res@trXMinF          = 1
 res@trXMaxF          = 4

 ;--set non-linear Gtw on left Y axis
 res@xyYStyle = "Irregular"
 res@xyYIrregularPoints = G_scale
 res@trYMinF          = Gtw_min
 res@trYMaxF          = Gtw_max

 ;resources for "right" Y axis
 res2 = True
 res2@tiYAxisString = "T~B~W~N~~S~day~N~ (~S~o~N~C)";+" [dashed]"
 res2@trYMinF         = Tw_min
 res2@trYMaxF         = Tw_max
 res2@xyLineColors      = "deepskyblue2"
 res2@xyDashPatterns    = 0; 3
 res2@xyLineThicknesses  = 1.5
 res2@xyLineOpacities    = 1
 res2@pmLegendDisplayMode    = "Never"

 res2@tmYRMode        = "Manual"	
 res2@tmYRTickStartF  = 26
 res2@tmYRTickEndF    = Tw_max; 38;
 res2@tmYRTickSpacingF= 3
 res2@tmYRMinorPerMajor=2
 res2@tmYRAutoPrecision=False
 res2@tmYRPrecision   = 2
 ;res2@tmBorderThicknessF    = 1.5

 res@tfPolyDrawOrder    = "Predraw"               ; put line on top
 res2@tfPolyDrawOrder   = "Predraw"               ; put line on top


 ;-----------Fig3----------
 res@vpYF          = 0.90   ; location of plot1
 res@vpXF          = 0.10    ; location of plot1
 res2@vpYF          = 0.90   ; location of plot1
 res2@vpXF          = 0.10    ; location of plot1

 ;=====make sure TWday and Gtw overlap======
 res@xyDashPatterns    =  res2@xyDashPatterns
 res@xyLineColors      =  res2@xyLineColors
 res2@gsnStringFontHeightF = 0.014; for size of (a)(b)

 res@xyDashPatterns        = 3;
 ;res2@gsnCenterString      = "Hottest 1% of land grid cells"
 res2@gsnLeftStringOrthogonalPosF = 0.03 ;adjust the space/distance of gsnLeftString/gsnCenterString
 res2@gsnCenterStringOrthogonalPosF = 0.03
 plot1 = gsn_csm_xy2(wks1, xx, Gtw99g_median, TWday99g_median, res, res2); left Y for Gtw, right for TWday

 gsres                   = True                        ; poly res
 gsres@gsFillOpacityF    = 0.1
 gsres@gsFillColor       = "deepskyblue"                 ; color chosen
 dummy2 = gsn_add_polygon (wks1,plot1@xy2,xp,yp4,gsres) ;add TWday99g to the scalar id "xy2"


 ;=============G under full radiation========
 res@xyLineColors          = "purple2";
 res@xyDashPatterns        = 0;
 res@xyLineThicknesses     = 1.5
 res@xyLineOpacities       = 1
 plot1_0 = gsn_csm_xy(wks1, xx, Gday99g_sun_median, res);

 gsres@gsFillColor       = "purple2";                ; color chosen
 dummy1 = gsn_add_polygon (wks1,plot1_0,xp,yp1,gsres) ;Gday99g_sun
 overlay(plot1,plot1_0)

 ;==============under only diffuse radiation==========
 res@xyLineColors        = "orange2"
 res@xyDashPatterns     = 0; 1
 plot1_1 = gsn_csm_xy(wks1, xx, Gday99g_diff_median, res); 
 gsres@gsFillColor       = "orange"                    ; color chosen
 dummy3 = gsn_add_polygon (wks1,plot1_1,xp,yp2,gsres) ;Gday99g_diff
 overlay(plot1,plot1_1)

 ;=====under zero solar radiation==========
 res@xyLineColors        = "slategray"
 res@pmLegendDisplayMode = "Never"
 plot1_2 = gsn_csm_xy(wks1, xx, Gday99g_zero_median, res);
 gsres@gsFillColor       = "slategray4"                    ; color chosen
 dummy4 = gsn_add_polygon (wks1,plot1_2,xp,yp3,gsres) ;Gday99g_zero
 overlay(plot1,plot1_2)

  ;--add dark scenario without sweat limit but with moderate MH=300/1.7
  res@xyDashPatterns     = 1
  plot1_3 = gsn_csm_xy(wks1, xx, Gday99g_zero_median_noSw, res);
  overlay(plot1,plot1_3)

  ;--add dark scenario without sweat limit and resting (MH=100/1.7) --should overlap with TWday99g
  res@xyDashPatterns     = 2
  plot1_4 = gsn_csm_xy(wks1, xx, Gday99g_zero_median_noSw_rest, res);
  overlay(plot1,plot1_4)

 

  ;=add 0 reference line
  yvalues = (/0,0/)
  res@gsnYRefLine      = yvalues
  res@gsnYRefLineColor = "grey"
  res@gsnYRefLineDashPattern = 1
  res@gsnYRefLineThickness   = 1
  res@pmLegendDisplayMode    = "Never"
  plot01 = gsn_csm_y(wks1,yvalues,res) ; use same res as base plot
  overlay(plot1,plot01); plot01 becomes part of plot1
 

  ;--add some text
  txres                     = True                 ; text mods desired
  txres@txFontHeightF       = 0.012                ; default size is HUGE!
  txres@txJust              = "CenterLeft"         ; puts text on top of bars
  ;gsn_text(wks1,plot1,"Hottest 1% of land grid cells", 1.2, res@trYMaxF*0.95,txres) ;

  ;=add 0 reference line
  yvalues = (/0,0/)
  res@gsnYRefLine      = yvalues
  res@gsnYRefLineColor = "grey"
  res@gsnYRefLineDashPattern = 1
  res@gsnYRefLineThickness   = 1
  res@pmLegendDisplayMode    = "Never"
  plot01 = gsn_csm_y(wks1,yvalues,res) ; use same res as base plot
  overlay(plot1,plot01); plot01 becomes part of plot1
 

  ;--add some text
  txres                     = True                 ; text mods desired
  txres@txFontHeightF       = 0.012                ; default size is HUGE!
  txres@txJust              = "CenterLeft"         ; puts text on top of bars
  ;gsn_text(wks1,plot1,"Hottest 1% of land grid cells", 1.2, res@trYMaxF*0.95,txres) ;

  ;--add legend for shade
  lbres                    = True
  lbres@vpWidthF           = 0.15          ; labelbar width
  lbres@vpHeightF          = 0.075;0.1          ; labelbar height
  lbres@lbBoxMajorExtentF  = 0.6           ; smaller value for more space between color boxes
  lbres@lbMonoFillPattern  = True
  lbres@lbPerimOn          = False
  lbres@lbBoxLinesOn       = False
  lbres@lbFillOpacityF     = 0.2; 0.1
  lbres@lbFillColors       = (/"slategray4","orange","purple"/) 
  labels = (/"","",""/); no labels for the shade bars
  gsn_labelbar_ndc(wks1,3,labels,0.098,0.68,lbres) ; draw left labelbar column

  ;legend for lines, more general function gsn_legend_ndc
  lgres                    = True
  lgres@vpWidthF           = 0.125                ; width of legend (NDC)
  lgres@vpHeightF          = 0.075;0.1           ; height of legend (NDC)
  lgres@lgLineThicknessF   = 1.5                 ; line thicknesses
  lgres@lgPerimOn          = False               ; no perimeter
  lgres@lgLineColors       = (/"slategray","orange2","purple2"/)
  lgres@lgDashIndexes      = (/0,0,0/)          ; line patterns
  labels := (/"G~S~day~N~ (dark)","G~S~day~N~ (shade)","G~S~day~N~ (sun)"/)
  lgres@lgLabelOffsetF     = 0.08
  lgres@lgLabelFontHeightF = 0.07;          ; font height (also scales with vpWidthF)
  gsn_legend_ndc(wks1,3,labels,0.14,0.68,lgres)


  lbres@lbFillColors       = (/"deepskyblue","white","white"/) 
  labels = (/"","",""/); no labels for the shade bars
  gsn_labelbar_ndc(wks1,3,labels,0.098+0.16,0.68,lbres) ; draw left labelbar column
  lgres@vpWidthF           = 0.13           ; width of legend (NDC)
  lgres@vpHeightF          = 0.075           ; height of legend (NDC)
  labels := (/"T~B~W~N~~S~day~N~ (G~B~TW~N~~S~day~N~)","G~S~day~N~ (dark**)", "G~S~day~N~ (dark*)"/)
  lgres@lgLineColors       = (/"deepskyblue2","slategray","slategray"/)
  lgres@lgDashIndexes      := (/0,2,1/)          ; line patterns
  gsn_legend_ndc(wks1,3,labels,0.14+0.16,0.68,lgres)

  draw(plot1)
  print("finish plot 3")




end 
